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I. ABSTRACT 

In this study we have incorporated two time scales into the phase field crystal model 
of a binary alloy to explore different solute trapping properties as a function of crystal- 
melt interface velocity. With only diffusive dynamics, we demonstrate that the segregation 
coefficient, i^' as a function of velocity for a binary alloy is consistent with the model of 
Kaplan and Aziz where K approaches unity in the limit of infinite velocity. However, with 
the introduction of wave like dynamics in both the density and concentration fields, the 
trapping follows the kinetics proposed by Sobolev |23], where complete trapping occurs at 
a finite velocity. 

II. INTRODUCTION 

Due to the relationship between complex dentritic microstructures and the mechanical 
properties of welded, soldered and cast components, solidification of alloys is one of the 
most researched subjects in materials science. Computational modeling of the solidification 
process has been an integral part of understanding the underlying physics that is associated 
with all aspects of solidification. In the past few decades, modeling of rapid solidification 
and solute trapping has gained attention for its role in many modern solidification processes 
such as thermal spraying, spin coating and laser melting [UES]. As the solid-liquid interface 
position advances at very fast rates (~ 1 m/s), there is substantial aberration from the 
equilibrium concentrations at the interface in both the solid and liquid phases. This deviation 
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can be characterized through the segregation coefficient K, which is defined as the ratio of 
the sohd concentration Cg to the hquid concentration Ci at the interface. The increase in 
the segregation coefficient with increasing interface velocity V is known as solute trapping 
and a complete understanding of the trapping process requires the dependence K{V) on the 
kinetic and thermodynamic properties of the alloy. 

There are three main theories of solute trapping. The most widely known is the continuous 
growth model (CGM) of Aziz and coworkers [21 |3]. The CGM introduces the so-called 
diffusive speed, Vd, which is the velocity at which a solute atom can traverse the solid- liquid 
boundary. Aziz and Kaplan define this parameter as Vd = D/X, where D is the diffusion 
coefficient in the interface region and A is the width of the interface. According to the CGM, 
significant trapping occurs when the interface velocity becomes comparable to Vd and the 
main prediction of the model in the limit of dilute concentrations is given by: 

-(-'^^ « 

where k^, is the equilibrium segregation coefficient. Notice, that the complete trapping limit 
given by K{V) = 1 is approached asymptotically as V — i- oo. The work of Jackson et 
al |16j demonstrates similar asymptotic behavior, but K{V) is described as a power law 
given by: 

K{V) = (2) 

Where A is a parameter analogous to the reciprocal of Vd- 

In contrast, the work of Sobolev ^31 predicts that there is an abrupt change in the 
segregation coefficient at a finite velocity. In other words, complete trapping {K{V) = 1) 
occurs at a well defined velocity. Their theory can be explained in terms of a "caging" 
phenomenon in the bulk liquid of the system. A solute atom spends a short amount of time 
(~ pico seconds) in a cage formed by nearest-neighbor atoms before it undergoes a random 
hop to another cage-like shell. If the interface velocity V is sufficiently fast, the liquid atoms 
attach themselves to the solid crystal at a rate that is comparable to this short time scale. 
At very high velocities the atoms do not have enough time to escape their cage and they 
completely freeze in their location. To capture the short time scale, Sobolev modified Pick's 
first law of diffusion by introducing the following: 

B T 

J + TB-^ = -DVC (3) 
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where tb is the relaxation time required for the flux to reach its steady-state regime. The 
fast atomic interactions in the liquid are incorporated into the time scale tb, which can be 
defined as tb = where is the bulk hquid diffusive speed. describes how fast 
a solute atom can travel in the liquid. Sobolev emphasizes that solute atoms have a finite 
velocity in the liquid, which differs from the Aziz's model where the diffusion speed in the 
bulk liquid is assumed infinite. Combining Eq. [s] with the conservation law ^ = —V ■ J, 
the diffusion equation in a steady state co-moving reference frame {x = Xo — Vt) can be 
written as: 

The solution of Eq. |4] for the solute concentration ahead of the interface (x > 0) is: 



CW = (a-C.)exp(-j^— |^^)+C„ V<Vg 

C{x) = C„ V> Vg (6) 

When X = (interface position), the concentration is defined as Cj. As x — )■ oo (far 
bulk liquid position), the far field concentration is Co- The solution given by Eq. |5] is 
similar to the well known solution of the concentration field in the liquid ahead of a moving 
interface except now the diffusion coefficient can be replaced by effective diffusion coefficient 
D* = D[l — (y/V^y]. For the detailed solution to Eq. |4| the reader should consult the 
following references [23l l2l] . 

The above model extends the Aziz model Eq. [T] by substituting the diffusion coefficient 
with the effective one, D*, that was obtained from Eq. |5| Therefore, the segregation co- 
efficient in this case reaches unity when V > . The Sobolev/Galenko model has the 
following formulation: 



KiV) = 1 V>VE (6) 



As the velocity of the interface increases signifcantly, the atoms in the liquid are not able 
to escape from the solid back to the liquid and so they freeze in their positions as the system 
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solidifies. In other word, the solute atoms do not have enough time to jump back to the 
liquid once the interface velocity equivalent to the bulk diffusive velocity, which results in 
complete trapping. 

There is mounting evidence that the Sobolev description of solute trapping is an accurate 
one. Danilov and Nestler [8] showed that the Sobolev prediction was a better fit for the 
experimental data on Si- As alloys by Kittl et al [18], especially in the high velocity limit. 
More recently, splat cooling experiments of Al-Mg were performed by Galenko and Her- 
lach [IS] and the results show a change from a eutectic to super saturated solid solution at 
a finite velocity. Furthermore, theoretical work of [I2j shows that the interface velocity as a 
function of undercooling changes from a power law to a linear relationship due to complete 
trapping and the model shows a very good agreement with Cu-Ni experiments. In addition 
to experimental studies, early Molecular Dynamics (MD) work by Cook and Clancy showed 
partitionless crystal growth (-^'(V") ^ 1) in a Leonard- Jones system at a finite velocity [7]. 
Recently, Yang et al also showed that complete trapping does occur at finite velocity. The 
authors simulated a Leonard- Jones binary and a Cu-Ni EAM model to show that the results 
hold for different systems [28]. It was found that the Sobolev model was a better fit to their 
MD data, while the Aziz model underestimated the high velocity cases when K{V) = 1. 
The Yang et al. results confirmed complete trapping but MD simulations are atomistic 
and cannot show whether there is a need for the second time derivative in the continuum 
description that was suggested by Sobolev [23] • 

The phase field models of Wheeler et al [T] and Echebarria el al [5] do not show com- 
plete trapping in the high velocity limit using the parabolic form of the dynamic equations. 
Galenko et al modified these models to their hyperbolic form i.e. adding the second time 
derivative to the concentration field [13]. They observed that K{V) does indeed tend to 
1, however its behavior varies depending on the phase field phenomenology they choose. 
In addition, the diffusion coefficient of these models is a function that can be tailored to 
produce different trapping behaviors. 

In this paper we study solute trapping using a phase field crystal (PFC) alloy model with 
inertial dynamics. The PFC model captures the atomic scale structure analogous to MD 
yet its diffusive time scale is of the same order as regular phase field models. Also, the PFC 
is motivated from classical density functional theory and contains only a minimal number 
of parameters, which in principle, can be derived from fundamental liquid state properties 



5 



through the direct correlation functions controlhng the excess energy. These features are 
relevant to solute trapping since most of the interactions occur on the liquid side of the 
solid-liquid interface. Moreover, changing the parameters and examining their effects on 
the system is rather fast compared to MD. For instance, the mobility of the system can 
be adjusted without the need for new interatomic potential(s) in order to see the effect of 
the change on the trapping behavior. The work of Stefanovic et al pOf [21] was the first 
to add a second time derivative to the pure PFC model [l7j. They motivated the second 
time derivative from hydrodynamics and used the resulting kinetic equations to study the 
deformation and plasticity in nano-crystals. The addition of the second time derivative to a 
two component alloy PFC model will be discussed in the model section of this paper. The 
competition between inertial dynamics and mobility in both the concentration and density 
field dictates the trapping behavior of the solute atoms will be examined. 



The rest of the paper is organized as follows. Section III describes the PFC model free 
energy as well and the equations of motion for the density and concentration. Section IV 
presents our findings by examining the relationship between the wave terms, the mobility of 
the system and their effect on the concentration and density fields. We end the paper with 
conclusions in Section IVl 

III. MODELING APPROACH 
A. PFC Model 

For a binary alloy system that consists of A and B atoms we start by introducing the 
PFC free energy functional jTT] : 
F 

The fields n and ip are the dimensionless local number density and concentration, respectively 



rfr{i?4 + (2V^ + V^) n - + + + + (7) 



and they can be approximated as follows: 

PA + PB ^ I PA-PB .gx 

n = 1 ip = z (8) 

Pe P 

Pa is the number density of species A (Solvent) and is the density of species B (Solute), 
p is the average density, pe is a reference liquid density, ks is Boltzmann's constant, R is the 
average atomic radius, which sets the length scale of the system, and T is the temperature. 
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The last three terms in Eq. [7] represent a free energy of the Cahn-Hilliard type [6], which has 
been used extensively to study phase separation. is the dimensionless bulk modulus of the 
liquid and is the dimensionless bulk modulus of the solid, which control the energy scale 
of the system. Following [TT] they are expanded as follows Bg = B^ + -82^'^ and B^ = Bq. 
The difference between the bulk moduli Ai?o = Bq — Bq sets the temperature scale of 
the system. z/,t and u are constant parameters related to the direct two point correlation 
functions and they can be calculated using classical density functional theory or can be fitted 
to phenomenological models to describe various materials properties. It is worth noting that 
minimization of the free energy yields hexagonal structures in 2D, BCC in 3D and a constant 
density profile in the liquid. The model is capable of producing two different phase diagrams, 
eutectic and spinodal binary systems. Eq. [7] is a reasonably simple model that can be used 
to simulate solidification, phase segregation and elasticity/plasticity [22] . 

Diffusion-controlled solidification dynamics is the PFC formalism are simulated via: 



V-(m,V^^\ (9) 



dt \ 5pB, 

where and Mb are the mobilities of each atomic species. The dynamics is driven by 
minimizing the free energy of the system in similar fashion to the Cahn-Hilliard equation [6] . 
One can rewrite Eqs. |9] in terms of the density n and the concentration fields using 
Eqs.gdl]: 



— = V ■ MiV— + V ■ M2V— 
ot on dip 

?l^^.M,¥-^ + V.M,v'-f (10) 

ot dip dn 

where Mi = (Ma + MB)/pj and M2 = (Ma - Ms) /pi- Notice, we will neglect the effects 
of thermal fluctuations in all of the simulations. For substitutional diffusion we assume the 
mobilities of species A and B are the same, which implies that M2 = and the cross terms 



vanish from Eqs. 10 



As mentioned earlier, Sobolev [23] introduces a second time derivative in concentration to 
capture the fast interactions in the bulk liquid. In the case of PFC alloys, Stefanovic et. al 
considered inertial dynamics in both the concentration and density field p6j. They derived 
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these as an extension of this concept in pure materials [20] to binary alloys. Equations of 
motion for the concentration and density in this limit are given by: 



where M = Mi, the /3 and 7 coefficients are responsible for introducing two time scales in 
the dynamics of n and tjj. The ratio of jj and jj be shown to be important quantities 
that will be discussed in depth in the next section. Each ratio controls the effect of atomic 
inertial interactions in its respective field. 



The simulation domain is a 2D channel consisting of 4096 grid points in the (x) direction 
normal to the solid-liquid interface and 64 grid points in the parallel direction (y). To 
initiate the simulation, a solid seed (periodic density) of (32a; x 6Ay was positioned at the 
far left part of the simulation box while the rest of the simulation box was filled with liquid 
(constant density). All simulations were performed with periodic boundary conditions in 
both the X and the y directions. The average concentration was set to be tpo = —0.175 and 
the temperature ABq = Bq — Bq = 0.08. As mentioned in the introduction, the Aziz model 
(Eq. [T| and the Sobolev (Eq. [6]) model are derived for the dilute limit. The parameters 
above are thus chosen to ensure the system has a dilute concentration of B atoms and the 
average concentration is such that it lies in the solid-liquid coexistence region. 

Simulations were done for a eutectic phase diagram shown in Fig[T] Since the model only 
supports hexagonal symmetry in 2D and the simulation box is of a rectangular shape, some 
stress is present in the system, which can cause a change in the equilibrium concentration. To 
account for this, and any numerical inaccuracies associated with analytically calculating the 
phase diagram, we numerically evaluated the coexistence lines that were used to conduct our 
simulations. Analytical and numerical coexistence lines of the PFC alloy model are shown 
in Fig. ([T|. In particular, the black lines represent the analytical solution obtained from an 
amplitude expansion technique as described in [H] while the red lines show the the equilib- 
rium solid/liquid concentrations obtained from the simulations. The equilibrium simulation 
was conducted using the simulation setup mentioned in the previous paragraph. Increasing 




(11) 



B. Numerical Methodology 



8 




FIG. 1: Eutectic phase diagram plot of Ai?o vs ipo for the parameters Bq = 1.00, B2 = —1.80, 
t = 0.60, = 1.00, u = 4.00 and uj = 0.008. The insert shows the difference between the analytical 
and simulated coexistence lines of the phase diagram. 

the system size in the y-direction relieves most of the stress in the system but a plateau is 
reached after roughly 64 grid points, with the simulation time increasing significantly for 
larger systems without a considerable amount of further stress relief. The remaining dis- 
crepancy between the analytic solution and the full numerical simulation shown in Fig [T] is 
most likely due to the approximations involved in the amplitude expansion technique. It is 
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noted that the concentration on the left hand side of the phase diagram has a negative value, 
which implies the absolute value of the solidus concentration to be larger than that of the 
liquidus. In the PFC definition of the concentration (Eq. [S]) a value of = corresponds 
to an actual concentration of 0.5. Therefore in the results presented below, we will add 1 to 
both the solid and liquid interface compositions when reporting the segregation coefficient 
i.e. K{V) = {i;s + l)l{i^i + l). 

The method of driving the solid-liquid interface employed in the current study involves 
the application of a linear temperature (ASq) gradient in the direction normal to the bound- 
ary. The temperature gradient (G) is chosen such that the liquid phase remains above the 
equilibrium liquidus temperature and the solid is maintained at temperatures below the 
liquidus. Over the course of the simulation, the gradient is translated in the x-direction at 
the specified velocity. Some care must be taken in applying the temperature gradient. The 
temperature gradient choice affects the solidification transient regime, either prolonging or 
shortening it, depending on the steepness of temperature gradient. Moreover, a tempera- 
ture cap, that is, the horizontal portions of the temperature profile shown in Fig [2| in the 
bulk must be incorporated. Certain temperature caps can be used to ensure that the tem- 
perature gradient does not "outrun" the solid-liquid interface. If the cap is too small, the 
temperature of the solid is too close to the starting temperature, which gives the interface 
limited driving force. Therefore, the temperature gradient will outrun the interface at large 
velocities, which causes the system to solidify isothermally. In order to pull the interface at 
very high velocities the interface and the solid portion of the system must be assigned lower 
temperature caps. Additionally, the temperature cap prevents entering the AB part of the 
phase diagram, which causes the formation of a and {3 phases (i.e. precipitates). For the 
purpose of this paper the temperature cap and the temperature gradient were chosen to be 
A5o = ±0.03 and G = 0.0002 respectively. 

IV. RESULTS AND DISCUSSION 

A. Inertial Dynamics in the Concentration Field 

We start the discussion of solute trapping by examining the steady-state concentration 
profiles. To produce a ID concentration profile, the concentrations are averaged in the y- 



10 




FIG. 2: A typical 2D density profile with a temperature gradient and cap overlaid on top to 
illustrate the deriving force mechanism employed in the simulations. 

direction parallel to the interface. As will be seen in the results presented below, the solute 
concentration profiles exhibit periodic oscillations in the solid phase on a scale of its inter- 
planar spacing. The periodicity in the solid part of the concentration profile is a result of the 
density-concentration coupling term in Eq. [7} The average solid concentration in the steady- 
state regime can be calculated by averaging over the atomic planes of the solidified crystal. 
Refer to the appendix |X] for further explanation on the concentration profile calculations. 

Fig|3^ and[3|D show typical concentration profiles for slow {V = 0.02) and fast {V = 0.5) 
velocities, respectively, using diffusive kinetics(/3 = 0, 7 = 0). The solute accumulation ahead 
of the interface decreases as the velocity increases. It is clear that even in extremely fast 
velocities (FigjsjD) the concentration ahead of the interface is higher than the solid concentra- 
tion. When the second time derivative is activated in the concentration field (7/M = 100.0) 
and /3 = 10, there is a substantial amount of trapping in the low velocity profile as illus- 
trated in Fig [3]:. One can define the occurrence of complete solute trapping when the liquid 
concentration ahead of the interface equals the averaged solid concentration through the in- 
terface. Fig[3]i displays a concentration profile where complete trapping occurs, it is evident 
the solute peak ahead of the interface vanishes. 

A summary of the PFC model parameters is given by Table |T] for all the simulations 
reported in this study. The values of /3 and 7 were chosen so the ratios of these quantities 
give insight in the solute trapping behavior. One set of simulations has diffusive dynamics 
in both fields to provide a baseline for the effect of the inertial dynamics on the trapping 
behavior. 

As the velocity of the interface increases using diffusive dynamics in the concentration 
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FIG. 3: Typical concentration profiles; a and b simulated with diffusive dynamics, while c and d 
simulated with inertial dynamics. Frames e and f illustrate the corresponding density profiles at 
(e) V = 0.2 and (f) V = 0.5, both simulated with diffusive dynamics. 

and density, K{V) approaches unity in an asymptotic manner analogous to the Aziz model 
(7 = 0, /3 = (black in color data )) in Fig |4] or Fig [s] For the same set of parameters, 
but the second time derivative is activated in both fields, complete solute trapping can be 
achieved for a steady state front propagating at a finite velocity, as seen by the 7/M = 100 
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/3 7 ^ M M i Comments 

0.00 0.00 1.00 0.00 0.00 - Diffusive dynamics in both fields 
10.00 110.0 11.0 0.91 10.0 11.0 Wave dynamics in both fields 

14.0 2000.0 20.0 0.71 100.0 142.9 Wave dynamics in both fields 

10.0 2000.0 20.0 0.50 100.0 200.0 Wave dynamics in both fields 

14.0 200.0 20.0 0.71 10.0 14.29 Wave dynamics in both fields 

10.0 200.0 20.0 0.50 10.0 20.00 Wave dynamics in both fields 

10.0 100.0 20.0 0.50 5.0 10.00 Wave dynamics in both fields 

10.0 66.66 20.0 0.50 3.33 6.66 Wave dynamics in both fields 

10.0 33.33 20.0 0.50 1.66 3.33 Wave dynamics in both fields 

10.0 16.66 20.0 0.50 0.83 1.66 Wave dynamics in both fields 

TABLE I: Summary of simulation parameters in Eqs. [TT] 

and (3/M = 0.5 data (red in colour) in Figs. |4]or [s] We also found analogous behaviour, 
leading to complete trapping, was exhibited with /3 = 0. Figure |4] also shows that, for fixed 

, the amount of trapping rapidly increases with increasing . These observations suggest 
that complete solute trapping can not take place in the alloy PFC model without the second 
time derivative being present in the concentration field. The term 7 in the PFC context 
thus plays the same role as the tb parameter in Gelenko's theory. 

Conversely, it was found that for a fixed value of increasing ^ leads a reduction of 
trapping in the ip-field over the same range of velocities as that discussed above. This is 
illustrated in Fig |5} The 7/M = 100 data of Fig. [s] also suggest that an increase in /3/M 
requires a higher velocity to achieve complete trapping in the ip field. Lowering 7 to 10 has 
the effect of pushing to even larger values the trapping velocity (i.e. increasing V^). Thus, 
for lower values of 7, we expect that the velocity needed to achieve complete trapping (as f3 
increases) increases for lower values of 7. To better understand the role of the density field 
on the solute partitioning, it is instructive to consider the role that the inertial term of the 
density field plays in the in the kinetics of the interface. This will be examined in the next 
section. 
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FIG. 4: K(V) vs. Velocity with constant /3 and variable 7. The case of purely diffusive dynamics 
is also shown for comparison. The right panel is zoom in of the left panel. 
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FIG. 5: K(V) vs. Velocity with two values of 7 and variable (3. The case of purely diffusive 
dynamics is also shown for comparison. The right panel is zoom in of the left panel. 



B. Inertial Dynamics in the Density Field 

The behavior of diffusive density dynamics in rapid sohdification was examined by [27] 
using the pure PFC model of Ref. [17]. They show that as the velocity of the interface 
increases, the difference between the solid and liquid densities decreases. They demonstrate 
that the density field follows an Aziz type trapping, analogous to that of the concentration. 
They attribute this change in the density to vacancy trapping as the interface advances in 
rapid solidification. We observe similar behavior in the density field of the binary alloy, where 
in the slow velocity regime, there is a noticeable jump in the density across the the solid- 
liquid interface (Figlsfe). As the the interface deviates from the local equilibrium condition 
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at higher interface velocity, the density jump significantly drops, as shown in Fig[3]F. 

The use of inertial dynamics in the pure material to examine the transition from a periodic 
to homogenous solution was first discussed by Galenko et al [H] . In their work, they inves- 
tigate the stability of the second order differential equation for both the Swift-Hohenberg 
and the PFC models. At low driving force (small interface speeds), the inertial terms can 
be neglected to retain a parabolic form of the dynamics, since the transition from one phase 
to another occurs on diffusive time scales. As the driving force increases, the inertial term 
in the dynamics (hyperbolic form) provides an extra degree of freedom to capture rapid 
kinetics. They compare the parabolic and hyperbolic solution for the same [high] driving 
regime and find the hyperbolic form of the equation leads to slower transition speeds. They 
ascribe the slowing down in the front speed to the increase in the relaxation time caused 



by the inertial kinetics. As the (3 term increases in Eq. [TT| the system takes a longer time 
to switch to diffusive dynamics. We also simulated isothermal quench simulations (i.e. no 
temperature gradient) and similarly found that the solid-liquid interface slows down as the 
effect of the inertial term increases, and when /3/M is of order 100 or higher, it was found 
that the interface virtually stops over the time of the simulation, leading to subsequent 
concentration pile up ahead of the interface, and lower solute trapping. 

The situation is different for the case of directional solidification used here for binary 
alloys, since conditions at the solid-liquid interface are driven by a temperature gradient that 
imposes the front velocity at steady state. In this case, for a given set of parameters in the 
free energy, as the pulling velocity increases the interface temperature is found to decrease 
i.e. slides down the temperature gradient towards its lower cap (see Fig. [2]), effectively 
increasing the driving force to keep the interface moving with the imposed higher speed. It 
was also found that a decrease in local interface temperature occurred for a fixed interface 
velocity when /3/M was increased. This is shown in Figure |6| which plots the interface 
temperature versus the pulling speed for two values of /3/M, showing that increasing either 
speed or (3 has the effect of decreasing the local interface temperature. 

To better understand the role of /3 on interface kinetics, we considered the lowest order 
sharp interface limit of the amplitude equation corresponding to the density equation used 
in this study. Details of the derivation are shown in Appendix |Bj This approach begins 
with PFC free energy coarse grained in terms of the amplitude of the density field, (p, where 
= in the liquid and = 0^ in the solid, where 0^ is given in Appendix [Bj This coarse 
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FIG. 6: The interface temperature vs. the velocity of the interface for two values of /3. 



grained free energy is given by Eq. Bl, and was first introduced in [TU]. The equation of 



motion of the density ampHtude, 0, is given by Eq. |B2[ for a co-moving reference frame. 
The addition of the inertial term in the dynamics studied here re-scales the coefficient of 



the highest gradient term, i.e. C in Eq. B2). In the Umit of small velocities, Appendix [B| 
derives a relation between the interface temperature (AS*) versus the steady state interface 
velocity V, reproduced here for convenience, 

90 ,0 V 



24 ZM^C 

2 



(12) 



where C = W{h)'^ — f3V'^/M, = dx, is the steady state solution of 

the phase field equation and W{n) the anisotropy interface width of (p. The precise form 
of (j) is not known, however we assume assume that it exhibits a transition from to 1 
across the solid-liquid interface, analogous to a hyperbolic tangent solution, which is the 



lowest order solution of Eq. |B2| in the V = limit. Numerical solutions of the amplitude 
equation are consistent with this assumption. As shown in appendix |B] the assumption of a 
(1— tanh(?/))/2 profile leads to = 1/3. The variables ipi, i/jg are the interface concentrations 
on the liquid/solid sides of the interface. 



Figure M plots AB^ predicted by Eq. 12 using the parameters from our simulations for 



two values of /3/M. Eq. 12 has qualitatively similar shape and order of magnitude in the 



low V range as the data from the numerical simulations. The discrepancy in ABo for 



large velocities is a consequence of the crude low-l^ approximations used to derive Eq. 12 
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FIG. 7: Interface temperature vs. interface velocity using Eq. B9 The il^i and ips in Eq- 12 



were 



extracted from the numerical simulations of the corresponding 7 and /?. (a^ = 1/3 in Eq. 12). 



Despite the approximate form of Eq. 12, it allows us to understand the basic role of (3 on 
interface kinetics. For example, it predicts that increasing P at large V leads to a decrease 
in dimensionless interface temperature. It is noted that while our simulations show that 
increasing (3 also increases -ipi (i.e. decreases k(V) in Figure. the magnitude of xjjf is too 



small to account for the change in ABl in Figs. Kjlor M (Eq. 12 ). 



Equation 12 is the \ow-V analogue of the sharp-interface relation derived by Aziz [2J for 
the continuous growth model. The intertial kinetics of the density equation are seen to 



essentially renormalize the dimensionless wave speed via C factor in the last term of Eq. 12 



It is noted that at vanishing velocities, C becomes essentially constant independent of V. As 
velocity increases, however, the role of /3 is enhanced. Specifically, Stefanovic et. al showed 
that /3 increases the effective time the system spends in wave like dynamics. This will have 



the effect of increasing the effective kinetic term in the sharp interface term in Eq. 12 



It is instructive to derive a more precise form of the interface kinetics for the alloy PFC 



model that couples information from the equation for ip with Eq. 12 This can be done 
more accurately than what was done here by solving the model equations perturbatively 
in a small variable like C, and matching solution in the interface regions to those outside 
the interface (i.e. bulk) regions. This will allow us to simultaneopusly relate the interface 
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temperature AB^ to both /3 and 7. This will be presented elsewhere. 

V. CONCLUSION 

The present study of rapid solidification using a binary PFC model was conducted to 
understand the role of the inertial dynamics in solute trapping in the PFC alloy formalism. 
According to Sobolev's prediction [231 [21], the hyperbolic form of the concentration diffusion 
equation predicts complete solute trapping to occur at a finite velocity. Our PFC alloy 
simulations are consistent with this theory, with the transition to complete trapping being 
controlled by the ratio (which controls inertial dynamics in the solute field). In particular, 
Increasing 7 reduces V^, leading to complete trapping at lower velocities. Since our alloy 
model couples the concentration field to that of the density, the role of rigidity in the density 
field was shown to modulate the role of inertial dynamics in the solute field. In particular, 
increasing the time scale over which the density exhibits inertial dynamics (by increasing the 
parameter /3 in the density equation), leads to a decrease in the interface temperature, and 
an increase in concentration on the liquid side of the interface, the latter of which decreases 
the partition coefficient. Our numerical results were consistent with a new sharp interface 
equation we derived to relate the interface temperature to the local interface velocity. 

Appendix A 

To compute the values of K{V) in Fig|5]and Fig|4| we need to calculate the averaged solid 
and liquid concentrations at the interface. However, in the solid phase, the concentration 
oscillates at a wave-length equivalent to the inter planar spacing. Therefore, to obtain the 
values of the averaged concentrations, we employ a Fourier filter to smooth the concentration 
profiles as shown below: 

fn = fne^-'^'^y' (Al) 

n is the grid spacing, /„ is the fourier transform of the concentration profile, /„ is the 
smoothed fourier transform and e controls the wave modes of the fourier filter. We then 
back transform the smoothed profile to real space /„. To calculate the minimum error 
associated with this method, we utilize a similar method used by [5l [19] where the error was 
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calculated by taking the second order derivative of the smoothed profile and sum all the 
terms over the length of the simulation box as shown in the following equation: 

N 

S = Y.^8Uf (A2) 

n=0 

where 5fn = (/n+i + /n-i — '^fn)/dx'^, which represents the second order derivative of the 
smoothed concentration profile. As the value of e increases, the error contributions from 
{5fnY decreases as the profile becomes smoother. 



We have modified Eq. A2 by adding error contributions from the interface liquid concen- 



tration in the following manner: 



N 



S = J2(^fn? + [Pa - Ps? (A3) 



n=0 

Pa and Ps are the liquid concentration peaks of the actual concentration profile and the 
smoothed smoothed profile, respectively. When e has a small value, the two peaks overlap 
on top of each, however, as the e value increases, the smoothed peak deviates from the actual 
concentration peak causing an increase in the 5* value. The two terms of Eq. |A3| have error 
contributions and they vary as we manipulate e. To minimize the error, we vary the value 
of e as illustrated in Fig|8] We found an optimal value of e = 0.343. 

The systematic error associated with the filter method (^ 6 x 10~^) is bigger than the 
statistical error resulted from the simulations. Since the noise terms were ignored for these 
simulations, there no is significant statistical error produced. After doing the error propa- 
gation analysis the final error associated with KiV) is ~ 1 x 10^^. 

Appendix B 

The derivation of the interface temperature begins with the free energy functional of our 
model, expressed in terms of an amplitude field, and concentration ip. The former is like 
an order parameter. The form of this free energy is derived by Elder et al jlUj : 

7- = I £xlw\n) ' + iB- SABA' - 4(/ + (27 - 
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FIG. 8: Error dependency on the value of e 

where W'^{n) is the interface width and all the other coefficients are related to PFC model 
and defined in section lllli 



The equation of motion for cj) corresponding to the density field in Eq. [TT| written in a 
steady-state co-moving reference frame, is given by: 



di]"^ M dr] d(j) 



(B2) 



where C = ^^^(n)-^ and ^ = 6[A5o + 5^^2]0- 12t02 + 9Oz/03^ the latter of which is the 



derivative of the non-gradient part of the integrand in Eq. Bl To simplify the analysis of 



Eq. B2 we will consider the low velocity regime, where we can assume to reasonable accuracy 
that the steady state solution of Eq. |B2| can be approximated by: 

V 



s9 



C 

<Ps dg{y) 



(B3) 



'Vc 



dr] ^JC dy 

where 0s = (t/15z/)(l + (a/1 — 15(Ai?o + B\iIjI)v /t"^)) is the far field value of in the solid 
[in], evaluated here at the capped value of Ai?o in the far field solid. The function g{ii) 
satisfies, giy — )■ ±oo) = 0, 1. While crude, this approximation allows us to elucidate some 
important features of the role of (3 in the interface kinetics of our model. 
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To extract the relationship between interface temperature and the interface velocity, we 



employ a projection method, whereby we multiply both sides of Eq. B2 by d(j)/dri and 
integrate from — oo to oo, which yields: 

I ' ( 

M 



oo J —oo 



df] 



'df 



dT] = 



(B4) 



drj drj^ ' ' M \dri J ' drj 

From the boundary conditions of g{y), it is straightforward to show that the first term 
vanishes. Using the second of Eq. |B3| the second term becomes 



M C 



dg{u) 



du 



du 



V 



(B5) 



where 



i^d/dy)'^ dy- example, if we employ g{y) 



tanh(y)) /2, 



1/3. Since the last term of Eq. B4 involves ABo{(f))), we approximate it by considering 



the limit where vC ~ ^ 1, which allows the lowest order form of amplitude gradient to 
be approximated by 

lim ^ ^ -(ps^iv) (B6) 
w<i dr] 



This is very approximate but allows the last integral of Eq. |B5| to be simplified as 

d(f) df 



dr] dcp 



drj 



r)=0 



using the definition above for df /d(f) and evaluating it at the interface, where 0(0) = ^ 
'0(0) = ipi (interface liquid concentration), results in 



6 + i?^0f] I - 12t (f)' + 90z. 



(B7) 
and 

(B8) 



Equations B5 and B8 are substituted into Eq. B4, yielding 



A5! 



90 



t<Ps - 5^0f - Y^^^ 



V 



(B9) 
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